本发明提供了一种基于局部方差分析的自适应曲波阈值地震去噪方法,该方法在深入分析曲波系数在不同尺度、不同角度下特征的基础上,提出了一种自适应阈值技术实现每个角度下有效信号与噪声的分离,并且借助曲波域的局部方差信息增强对非噪声成分的保护,从而减小对有效信号(特别是弱信号)的伤害。在合成地震记录与陆上地震数据的实践证明,阈值迭代法的实际去噪表现并不理想,而本发明去噪效果更加真实有效,在压制噪声与保护细节方面具有独特优势。
1.一种基于局部方差分析的自适应曲波阈值地震去噪方法,其特征在于,所述方法包括以下步骤:步骤1:对地震数据进行曲波变换获得曲波变换系数;步骤2:利用最精细尺度下曲波系数的标准差来确定随机噪声在时间域的标准差,采用Monta-Carlo分析获得曲波域各尺度及角度的变换权重,并以此计算曲波域中每一尺度及角度下的噪声标准差;步骤3:利用含噪信号在每一尺度及角度下的噪声标准差,结合局部方差分析,推算有效信号的标准差;所述的有效信号的标准差的计算公式为: σ s ( j , l ) = max ( σ y 2 ( j , l ) - σ a 2 ( j , l ) , 0 ) ]]>式中ntent="drawing" img-format="tif" inline="yes" orientation="portrait" wi="78"/>代表含噪信号的总体方差,j,l分别代表尺度与角度,σ_y(j,l)代表输入含噪信号的总体标准差,σ_a(j,l)代表噪声标准差;步骤4:利用上述获得的噪声与有效信号的标准差建立阈值,进行软阈值处理;步骤5:进行曲波逆变换,获得时间域的去噪结果。
2.根据权利要求1所述的基于局部方差分析的自适应曲波阈值地震去噪方法,其特征在于,步骤1中所述的曲波变换采用USFFT或Wrapping方法。
3.根据权利要求1所述的基于局部方差分析的自适应曲波阈值地震去噪方法,其特征在于,在估算有效信号在曲波域的标准差σ_y时,引入一种局部方差处理方法,对指定尺度与角度下的曲波系数矩阵,逐点计算含噪信号的总体方差,即ntent="drawing" img-format="tif" inline="no" orientation="portrait" wi="700"/>上式表示计算以目标点(m,n)为中心的一个局部窗口W内的方差作为该目标点的方差,进而利用权利要求1所述的公式能产生更准确的有效信号在曲波域的标准差。
4.根据权利要求1所述的基于局部方差分析的自适应曲波阈值地震去噪方法,其特征在于,步骤4中所述的阈值设为: T ( j , l ) = μ σ ^ a 2 ( j , l ) σ ^ s ( j , l ) ]]>其中j与l分别代表尺度和角度,μ为调节因子,ntent="drawing" img-format="tif" inline="yes" orientation="portrait" wi="179"/>与ntent="drawing" img-format="tif" inline="yes" orientation="portrait" wi="170"/>分别代表随机噪声和有效信号在不同尺度和角度下的标准差。技术领域
本发明属于石油勘探领域,特别涉及一种基于局部方差分析的自适应曲波阈值地震去噪方法。
背景技术
地震资料去噪工作贯穿于地震数据采集、处理和解释的全过程之中,而高信噪比则是地震资料数字处理的首要任务。随着地震勘探技术的发展,常规地震数据去噪手段已经越来越不能适应当前高精度勘探的需求,传统去噪技术的局限性促使地震数据去噪技术不断发展,去噪技术的内容也愈加丰富。曲波变换是一种具有方向性的多尺度变换,它在小波变换的基础上引入了方向性参数,能够对图像边缘部分的直线或曲线提供更加稀疏的表示。第一代曲波变换是在脊波变换的基础上引入函数局部化和频带剖析的思想发展起来的,能进一步表示多维信号中的更为普遍的曲线奇异性(Candès和xa0Donoho,xa02003)。由于其数字实现的大量数据冗余,Candes又提出了第二代曲波变换及快速算法(Candèsxa0等,xa02006a)。曲波变换起源于图像处理,与地震数据处理的结合则始于2004年SEG年会(Herrmann和xa0Verschuur,xa02004),xa0随着快速离散曲波变换的应用,使得曲波变换在地震数据处理中得到了迅速发展,目前在地震数据降噪与插值(Hennenfentxa0和Herrmann,xa02006,xa02007;xa0Kumar和Herrmann,xa02009;Wangxa0等,xa02010;刘国昌xa0等,xa02011),xa0多次波分离与面波压制(Herrmannxa0等,xa02007b;xa0Zhengxa0等,xa02010),波场外推(Linxa0和xa0Herrmann,xa02007)与地震偏移成像(Herrmannxa0等,xa02007a)等方面得到了愈发广泛的应用。曲波变换是为了解决小波变换识别线性奇异性边缘的不足而提出来的,具有各向异性和多尺度的性质,并且曲波变换基函数的支撑空间是狭长形的,相对于二维小波变换而言对曲线型边界可以进行更稀疏的表示。相比于傅里叶变换与小波变换,曲波变换能够取得绝佳的稀疏表示(Herrmann等,xa02007b)。稀疏性意味着曲波域中少量的大幅值曲波系数能够恢复大部分的有效信号,通常情况下利用小于总量5%的最大曲波系数就已经实现输入了数据主要轮廓的重建。与有效地震信号相比,随机噪声缺乏方向特性,通常体现为较小的曲波系数,因此利用有效地震信号与噪声在曲波系数的这种差异,在曲波域通过设定阈值压制(或去除)弱小的曲波系数,保留大曲波系数,然后执行逆曲波变换就实现了去噪过程,这是曲波阈值法去噪的基本出发点。目前,曲波阈值去噪技术中以阈值迭代法为代表的单一阈值方法忽视了曲波系数的特点,通常会使曲波系数过分稀疏化,不利于识别和恢复噪声背景下所隐藏微弱有效信号成分(张军华xa0等,xa02010)。阈值迭代法是曲波域去噪的一种经典手段,由Daubechies等人(2004)提出。设记录的(二维)地震数据的向量表示为ntent="drawing" img-format="jpg" inline="no" orientation="portrait" wi="62"/>xa0xa0xa0xa0xa0xa0xa0xa0xa0xa0xa0xa0xa0xa0xa0xa0xa0xa0xa0xa0xa0xa0xa0xa0xa0xa0xa0xa0xa0xa0xa0xa0xa0xa0xa0(1)其中,ntent="drawing" img-format="jpg" inline="no" orientation="portrait" wi="13"/>与ntent="drawing" img-format="jpg" inline="no" orientation="portrait" wi="14"/>为分别为有效地震信号与噪声。在此,用ntent="drawing" img-format="jpg" inline="no" orientation="portrait" wi="18"/>代表曲波变换矩阵,ntent="drawing" img-format="jpg" inline="no" orientation="portrait" wi="20"/>代表有效信号对应的曲波系数向量按(从大到小排列),则上式等同于ntent="drawing" img-format="jpg" inline="no" orientation="portrait" wi="81"/>xa0xa0xa0xa0xa0xa0xa0xa0xa0xa0xa0xa0xa0xa0xa0xa0xa0xa0xa0xa0xa0xa0xa0xa0xa0xa0xa0xa0xa0xa0xa0xa0xa0xa0(2)地震去噪问题等价于在曲波变换域中获得对ntent="drawing" img-format="jpg" inline="no" orientation="portrait" wi="20"/>最佳估计ntent="drawing" img-format="jpg" inline="no" orientation="portrait" wi="16"/>,可以表述为如下非线性优化问题ntent="drawing" img-format="jpg" inline="no" orientation="portrait" wi="206"/>xa0xa0xa0xa0xa0xa0xa0xa0xa0xa0xa0xa0xa0xa0xa0xa0xa0xa0xa0xa0xa0xa0xa0xa0(3)求解上述优化问题就是要反演一组1-范数最小的曲波系数,若该组系数满足经逆曲波变换之后与原始数据差的2-范数小于某一噪声参量,则作为最终输出。然而上述ntent="drawing" img-format="jpg" inline="no" orientation="portrait" wi="17"/>需要通过阈值降速的迭代过程实现,实际执行过程较慢。Elad等人(2005)年提出了相应的替代简化方程ntent="drawing" img-format="jpg" inline="no" orientation="portrait" wi="185"/>xa0xa0xa0xa0xa0xa0xa0xa0xa0xa0xa0xa0xa0xa0xa0xa0xa0xa0xa0xa0xa0xa0xa0xa0xa0xa0xa0xa0(4)由此,问题转变为ntent="drawing" img-format="jpg" inline="no" orientation="portrait" wi="16"/>的优化问题,即寻找使得ntent="drawing" img-format="jpg" inline="no" orientation="portrait" wi="89"/>成立的最大ntent="drawing" img-format="jpg" inline="no" orientation="portrait" wi="16"/>,通常利用如下的阈值迭代技术实现求解。ntent="drawing" img-format="jpg" inline="no" orientation="portrait" wi="200"/>xa0xa0xa0xa0xa0xa0xa0xa0xa0xa0xa0xa0xa0xa0xa0xa0xa0xa0xa0xa0xa0xa0xa0xa0xa0xa0(5)式中,ntent="drawing" img-format="jpg" inline="no" orientation="portrait" wi="20"/>表示软阈值函数,ntent="drawing" img-format="jpg" inline="no" orientation="portrait" wi="192"/>。计算过程中,每一次循环都伴随着阈值的改变,以此不断更新曲波系数,直到满足迭代终止条件。事实上,阈值迭代法把所有尺度及角度下的曲波系数作为一个长向量统一考虑,通过迭代手段确定合适的阈值之后进行截断,并对保留的大曲波系数进行逆变换就实现了去噪过程。这种手段忽视了不同尺度下曲波系数的不同特点,用同一个阈值进行处理就难免会丢失有效信息。
发明内容
本发明提供了一种基于局部方差分析的自适应曲波阈值地震去噪方法,在深入分析曲波系数在不同尺度、不同角度下特征的基础上,以针对性阈值为基础借助局部方差信息增强对曲波域非噪声信号方差的计算结果,并针对不同尺度与角度建立相适应的阈值,在恢复非噪声信号成分(特别是弱信号成分)具有明显优势。在数值模型测试与陆上地震资料的去噪应用实际表明,基于自适应阈值法的曲波域去噪能力更加真实可信,有效保留了有效信号特征,服务于地震资料后续处理与解释工作。本发明的上述目的是通过以下技术方案实现的:一种基于局部方差分析的自适应曲波阈值地震去噪方法,所述方法包括以下步骤:步骤1:对地震数据进行曲波变换获得曲波变换系数;步骤2:利用最精细尺度下曲波系数的标准差来确定随机噪声在时间域的标准差,并以此计算曲波域中每一尺度及角度下的噪声标准差;步骤xa03:结合含噪信号在每一尺度及角度下的标准差,分析局部方差,推算有效信号的标准差;步骤xa04:利用上述获得的噪声与有效信号的标准差建立阈值,进行软阈值处理;步骤xa05:进行曲波反变换,获得时间域的去噪结果。进一步的,步骤1中所述的曲波变换采用USFFT(即Unquallyxa0Spacedxa0Fastxa0Fourierxa0Transform)或Wrapping方法。本发明所提供的基于局部方差分析的自适应曲波阈值地震去噪方法的基本原理如下:为解决目前单一曲波阈值法的不足和缺陷,及增强保护有效信号的需要上,本发明借助局部方差信息,建立一种新型的自适应阈值去噪方法,其中阈值函数采用软阈值法,而阈值则设为ntent="drawing" img-format="jpg" inline="no" orientation="portrait" wi="125"/>xa0xa0xa0xa0xa0xa0xa0xa0xa0xa0xa0xa0xa0xa0xa0xa0xa0xa0xa0xa0xa0xa0xa0xa0xa0xa0xa0xa0xa0xa0xa0xa0xa0xa0(6)其中ntent="drawing" img-format="jpg" inline="no" orientation="portrait" wi="14"/>与ntent="drawing" img-format="jpg" inline="no" orientation="portrait" wi="9"/>分别代表尺度和角度,ntent="drawing" img-format="jpg" inline="no" orientation="portrait" wi="17"/>为调节因子,ntent="drawing" img-format="jpg" inline="no" orientation="portrait" wi="53"/>与ntent="drawing" img-format="jpg" inline="no" orientation="portrait" wi="53"/>分别代表(曲波变换域中)随机噪声和有效信号在不同尺度和角度下的标准差。随机噪声在时间域的标准差可利用最精细尺度下的标准差近似,进而可利用Monta-Carlo分析推算曲波域中每一尺度及角度下的噪声标准差(王海松和王伟,xa02009),而有效信号的标准差通常利用公式(7)确定ntent="drawing" img-format="jpg" inline="no" orientation="portrait" wi="248"/>xa0xa0xa0xa0xa0xa0xa0xa0xa0xa0xa0xa0xa0xa0xa0xa0xa0xa0xa0xa0xa0xa0xa0xa0xa0(7)式中ntent="drawing" img-format="jpg" inline="no" orientation="portrait" wi="22"/>代表含噪信号的总体方差。为突出增强阈值保留有效信号的能力,本发明在估算有效信号在曲波域的标准差ntent="drawing" img-format="jpg" inline="no" orientation="portrait" wi="22"/>时,引入一种局部方差处理方法:对指定尺度与角度下的曲波系数矩阵,逐点计算含噪信号的总体方差,即ntent="drawing" img-format="jpg" inline="no" orientation="portrait" wi="253"/>xa0xa0xa0xa0xa0xa0xa0xa0xa0xa0xa0xa0xa0xa0xa0xa0xa0xa0xa0xa0xa0xa0xa0xa0(8)上式表示计算以目标点(m,xa0n)为中心的一个局部窗口ntent="drawing" img-format="jpg" inline="no" orientation="portrait" wi="20"/>内的方差作为该目标点的方差,进而利用公式(7)能产生更准确的有效信号在曲波域的标准差。以下结合实施例与附图说明本发明具体实施方式。
附图说明
图1xa0是合成地震记录模型(a)及其含噪(20%)模型(b);图2是传统阈值迭代法对含噪记录模型的去噪结果(a)及剔除的噪声(b);图3是本发明(自适应阈值法)对含噪记录模型的去噪结果(a)及剔除的噪声(b);图4是某原始陆上叠后含噪地震剖面;图5是阈值迭代法对该陆上地震资料去噪后的剖面(a)及剔除的噪声(b);图6是本发明(自适应阈值法)对陆上地震资料去噪后的剖面(a)及剔除的噪声(b);图7是本发明(自适应阈值法)(a)及传统方法(阈值迭代法)(b)在A区域的去噪结果对比;图xa08是本发明(自适应阈值法)(a)及传统方法(阈值迭代法)(b)在B区域的去噪结果对比;图9是本发明所述的基于局部方差分析的自适应曲波阈值地震去噪方法的流程图。本发明与现有技术相比的有益效果为:1、本发明所述的基于局部方差分析的自适应曲波阈值地震去噪方法,借助局部方差信息增强对曲波域非噪声信号方差的计算结果,并针对不同尺度与角度建立相适应的阈值,避免了目前主流曲波阈值法造成的曲波系数过分稀疏化的特点,在识别和恢复非噪声信号成分(特别是弱信号成分)方面具有独特优势;2、本发明所述的基于局部方差分析的自适应曲波阈值地震去噪方法,在数值模型测试与陆上地震资料的去噪应用实际表明,基于自适应阈值法的曲波域去噪能力更加真实可信,有效保留了有效信号特征,服务于地震资料后续处理与解释工作。
具体实施方式
实施例1一种基于局部方差分析的自适应曲波阈值地震去噪方法,包括以下步骤:步骤1:对地震数据进行曲波变换获得曲波变换系数;步骤2:利用最精细尺度下曲波系数的标准差来确定随机噪声在时间域的标准差,并以此计算曲波域中每一尺度及角度下的噪声标准差;步骤xa03:利用含噪信号在每一尺度及角度下的噪声标准差,结合局部方差分析,推算有效信号的标准差;步骤xa04:利用上述获得的噪声与有效信号的标准差建立阈值,进行软阈值处理;步骤xa05:进行曲波逆变换,获得时间域的去噪结果。进一步的,步骤1中所述的曲波变换采用USFFT(即Unquallyxa0Spacedxa0Fastxa0Fourierxa0Transform)或Wrapping方法。进一步的,步骤3中所述的有效信号的标准差的计算公式为:ntent="drawing" img-format="jpg" inline="no" orientation="portrait" wi="248"/>式中ntent="drawing" img-format="jpg" inline="no" orientation="portrait" wi="22"/>代表含噪信号的总体方差,xa0j,lxa0分别代表尺度与角度,ntent="drawing" img-format="jpg" inline="no" orientation="portrait" wi="50"/>代表输入含噪信号的总体标准差,ntent="drawing" img-format="jpg" inline="no" orientation="portrait" wi="50"/>代表噪声标准差。进一步的,在估算有效信号在曲波域的标准差ntent="drawing" img-format="jpg" inline="no" orientation="portrait" wi="22"/>时,引入一种局部方差处理方法,对指定尺度与角度下的曲波系数矩阵,逐点计算含噪信号的总体方差,即ntent="drawing" img-format="jpg" inline="no" orientation="portrait" wi="253"/>上式表示计算以目标点(m,n)为中心的一个局部窗口ntent="drawing" img-format="jpg" inline="no" orientation="portrait" wi="20"/>内的方差作为该目标点的方差,进而利用权利要求3所述的公式能产生更准确的有效信号在曲波域的标准差。进一步的,步骤4中所述的阈值设为:ntent="drawing" img-format="jpg" inline="no" orientation="portrait" wi="125"/>其中ntent="drawing" img-format="jpg" inline="no" orientation="portrait" wi="14"/>与ntent="drawing" img-format="jpg" inline="no" orientation="portrait" wi="9"/>分别代表尺度和角度,ntent="drawing" img-format="jpg" inline="no" orientation="portrait" wi="17"/>为调节因子,ntent="drawing" img-format="jpg" inline="no" orientation="portrait" wi="53"/>与ntent="drawing" img-format="jpg" inline="no" orientation="portrait" wi="53"/>分别代表随机噪声和有效信号在不同尺度和角度下的标准差。图1(a)xa0为地震资料处理领域中一个著名的合成地震记录模型,包含背斜、向斜、角度不整合多种地震现象,在此基础上加入20%的白噪声(即最大值的20%作为白噪声的标准差),作为含噪模型如图1(b)所示。利用公式(5)所示的阈值迭代法对其进行去噪,经过10次迭代去噪结果与剔除的噪声分别如图2(a)与2(b)所示。根据信噪比计算公式:ntent="drawing" img-format="jpg" inline="no" orientation="portrait" wi="174"/>,其中ntent="drawing" img-format="jpg" inline="no" orientation="portrait" wi="18"/>与ntent="drawing" img-format="jpg" inline="no" orientation="portrait" wi="17"/>分别代表无噪与去噪后数据,阈值迭代法去噪后信噪比达到14.89。不过,从图2(b)中可以看出阈值迭代法明显地损失了部分有效信号。对图1(b)所示的含噪地震记录模型采用自适应阈值法去噪,结果如图3所示。与阈值迭代法相比,该方法提出的噪声成分更纯粹,明显降低了对非噪声成分的伤害,去噪后信噪比达到15.6。这种优势一方面根源于在曲波域按不同尺度、不同角度分别建立阈值处理,另一方面局部方差信息的使用则有助于识别非噪声成分。以我国东部地区某二维测线地震资料进行叠后去噪处理,原始地震剖面如图4所示。除基底下的大量噪声之外,噪声并非均匀分布在整个地震剖面上,而是呈现“团块状”分布。分别运用阈值迭代法与本发明(自适应阈值法)进行去噪,结果见图5与图6。受迭代算法中白噪声的前提假设,图5(b)中阈值迭代法剔除的“噪声”显得更均匀,统计特征也更符合白噪声的特点,但这与实际噪声资料情况不符。相比之下,本发明(自适应阈值法)产生的图6(b)中剔除的噪声则更具针对性,如地震剖面上信噪比较低的基底、岩石破碎带(A、B两区)分布的大量噪声得到有效移除,而高信噪比区域的噪声剔除量相对少,因而使有效地震信号特征得以保留。进一步对比噪声发育的A区与B区,如图7与8所示,不难发现阈值迭代法的去噪能力欠佳,岩石破碎带依然残留明显噪声,相比之下本发明(自适应阈值法)提出的去噪方法使得波形特征变化更平滑,能量突变点显著减少,具备更突出、更具针对性的去噪能力。